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ABSTRACT 

Context. In order to study the terrestrial-type planet formation during the post oligarchic growth, the initial distributions of planetary 
embryos and planetesimals used in N-body simulations play an important role. In general, most of these studies typically use ad hoc 
initial distributions based primarily on theoretical and numerical studies. 

Aims. We analyze the formation of planetary systems without gas giants around solar-type stars focusing on the sensitivity of the 
results to the particular initial distributions used for planetesimals and planetary embryos at the end of the gas phase of the protoplan¬ 
etary disk. The formation process of terrestrial planets in the habitable zone (HZ) and the study of their final water contents are also 
topics of special interest in this work. 

Methods. We developed two different sets of N-body simulations starting with the same protoplanetary disk. The first set assumes 
typical ad hoc initial distributions for embryos and planetesimals and the second set obtains these initial distributions from the results 
of a semi-analytical model which simulates the evolution of the protoplanetary disk during the gaseous phase. 

Results. Both sets of simulations form planets within the HZ. Using ad hoc initial conditions the masses of the planets that remain in 
the HZ range from 0.66M^ to 2.21 M^. Using more realistic initial conditions obtained from a semi-analytical model, we found that 
the masses of the planets range from LlSA/jg to 2.21Af®. Both sets of simulations form planets in the HZ with water contents ranging 
between 4.5% and 39.48% by mass. Those planets that have the highest water contents with respect to those with the lowest water 
contents, present differences regarding the sources of water supply. 

Conclusions. Comparing both sets of simulations we suggest that the number of planets that remain in the HZ is not sensitive to the 
particular initial distribution of embryos and planetesimals and thus, the results are globally similar between them. However, the main 
differences observed between both sets are associated to the accretion history of the planets in the HZ. These discrepancies have a 
direct impact in the accretion of water-rich material and in the physical characteristics of the resulting planets. 
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1. Introduction 


The formation process of terrestrial-type planets goes through 
different stages. In the classical core accretion model 
IILissauer & Stevensonll2007h . the first step in the planet forma¬ 
tion process is the sedimentation to the mid-plane of the dust 
of the disk and the formation of planetesimals. The mecha¬ 
nisms which lead to the formation of planetesimals are still un¬ 
der debate. However, in the last few years, three main scenarios 
for the formation of planetesimals have emerged: formation via 
coagulation-fragmentation cycles and mass transfer from small 
to large aggregates, growth of fluffy particles and subsequent 
compaction by self-gravity, and concentration of small particles 
(often called pebbles) in the turbulent g as, and gravitational fr ag- 
mentation of over dense filaments (see iJohansen et al.ll20T4l for 
a detailed discussion and references therein). After the forma¬ 
tion of the first population of planetesimals, these bodies grow 
quickly by mutual accretion until some of them (often called 
embryos) start to separate from the population of planetesi¬ 
mals. This process is known as planetesimal runaway growth. 
When embryos begin to gravitationally dominate the dynamic 
of the population of planetesimals, the runaway growth regime 
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switches to the oligarchic growth regime. In this regime, the 
gravitational excitation of the embryos over the planetesimals 
limits their growth, and embryos are the only population that 
contin ue growing by accretion of planetesimals (iKokubo & Idal 
Il998h . 

The interaction of the embryos with the gaseous disk could 
lead to the migration of these bodies along the disk, a phe¬ 
nomenon known as type I migration. In idealized isothermal 
disks type I migration is alw ays inward and at high migra¬ 
tion rates dTanaka e t al.ll200 2l). However, if more realistic disks 
are considered (iPaardekooper et al.l 1201011201 ll) type I migra¬ 
tion could be outward. This outward migration can also occur 
when the planet achieves a significant mass (~ lOMm) and cor- 
rotation torques become not negligible (iMasset et al.ll200^ . But 
even in isothermal disks, type I migration could be outward if 
full magnetohydrodyna mic disk (MHD ) turbulence is consid¬ 
ered (|Guilet et al.ll2013l) . More recentlv, ]Benftez-Llambav et al.l 
(I2OI5I) found that if the released energy by the planet due to 
accretion of solid material is taken into account this phenom¬ 
ena generates a heating torque which could significantly slow 
down, cancel, or even reverse, inward type I migration for plan¬ 
ets with masses < SM,^. Embryos continue growing by accretion 
of planetesimals and other embryos when th e distance between 
them becomes less than 3-4 mutual Hill radii (lAgnor & Asnhau^ 
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I 2 OO 4 I) due to the fact that the damping generated by the gas re¬ 
mains the embryos in nearly circular and coplanar orbits. 

It is worth mentioning that in the last years a new alter¬ 
native model was proposed for the formation of giant planet 
cores and terrestrial planets. This model proposes that the core 
of giant planets and terrestrial planets could be formed as seeds 
that accrete small particles, often called peb bles. lOrmel & Klahil 
(l2010l) . [Lambrechts & Johansenl(l2012h . and iGuillot et al.l ( 2014h 
showed that pebbles, which are strong coupled to the gas, could 
be accreted very efficiently to form planetary bodies. 

More recently, iJohansen et al.l (1201 5h showed that other 
structures of our Solar System, like the asteroid belt and the 
transneptunian belt, could also be explained by the pebble ac¬ 
cretion mechanism. 

Once the gas has dissipated from the system, and maybe 
some giant planet have formed, the terrestrial planet region con¬ 
tains a great number of planetary embryos which still evolve em¬ 
bedded in a remnant swarm of planetesimals. Embryos continue 
growing by accreting planetesimals and by stochastic collisions 
with other embryos in a regime known as pos-oligarchic growth. 
Finally, the terrestrial-type planets form and settle onto stable 
orbits on time scales of the order of 10^ years. 

In general terms, the numerical simulations that analyze the 
late stage planetary accretion after the gas disk dissipation re¬ 
quire physical and orbital initial conditions for planetesimals 
and planetary embryos, as well as for the gas giants of the sys¬ 
tem. These initial conditions are tipically selected ad hoc by 
man y works that invest i gate the process of planetary forma¬ 
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The study of planetary systems without gas giants re¬ 
sults to_b£_very_jr^resting.In_^ck_seTCral_obseryational stud¬ 
ies dCumming et al.l 2008: iMay or et al. 201 ih and theoretical 
works dMordasini et al. 2009t Miguel et alJ 201 ih haye sug¬ 
gested that the planetary systems consisting only of rocky plan¬ 
ets wo uld seem to be the most common in the Uniyerse. Re¬ 
cently, iRonco & de Elfal (l2014ll haye analyzed the terrestrial- 
type planet formation and water deliyery in systems without gas 
giants around Sun-like stars. In particular, the authors focused on 
low-mass protoplanetary d isks for a wide range of surface den¬ 
sity profiles. In that work. [Ronco & de Elfal (l2014l) assumed ad 
hoc initial conditions for the distributions of planetesimals and 
planetary embryos after the gas dissipation. 

Here, we present results of N-body simulations aimed at an¬ 
alyzing the process of planetary formation around Sun-like stars 
in absence of gas giants. In particular, the present study focuses 
on the sensitiyity of the results to the particular initial distribu¬ 
tion adopted for planetesimals and planetary embryos after the 
gas dissipation. To do this, we decided to carry out two different 
sets of numerical simulations starting with the same initial pro¬ 
toplanetary disk and modifying the particular initial distribution 
for planetesimals and planetary embryos. The first set of numer¬ 
ical simulations assumes ad hoc initial condition s for planetesi¬ 
mals and embryos like in iRonco & de Elf^ d2014l) while the sec¬ 
ond set of numerical simulations obtains initial conditions from 
a semi-analytical model, which is able to analyze the eyolution 
of a system of planetesimals and embryos during the gaseous 
phase. A comparatiye analysis between these two sets of numeri¬ 
cal simulations will allow us to enrich our knowledge concerning 
the process of planetary formation. Moreoyer, we consider that 
the present work will allow us to clarify our understanding about 
a correct selection of initial conditions for numerical simulations 
of formation and eyolution of planetary systems. 


This paper is therefore structured as follows. In Sect. 2, we 
present the main properties of the protoplanetary disk used in our 
simulations. Then, we outline our choice of initial conditions in 
Sect. 3 and discuss the main characteristics of the N-body code 
in Sect. 4. In Sect. 5, we present the results of all our simulations. 
We discuss such results within the framework of current knowl¬ 
edge of planetary systems in Sect. 6 and present our conclusions 
in Sect. 7. 


2. Protoplanetary Disk: Properties 

The properties of the prot oplanetary disk considere d in this work 
are the same presented in iRonco & de EH^ (l2014ll . The gas sur¬ 
face density profile that represents the structure of the protoplan¬ 
etary disk is giyen by 





( 1 ) 


where is a normalization constant, a characteristic radius 
and y the exponent that represents the surface density gradient. 

can be found by integrating Eq. [1] oyer the total area of the 
disk and is a function of y, and the mass of the disk Mj. 

Analogously, the solid surface density profile 2s(r) is repre¬ 
sented by 


2s(r) = S"?7ice 



( 2 ) 


where rjice represents an increase of a factor of 4 in the amount 
of solid material due to the condens ation of water beyond the 
snow line, which is located at 2.7 au (lHayashilll98lh . Although 
other authors fo und a much lower yalue for this increase of solids 
( LoddersI l2003h we used the same factor as IRonco & de Elfal 
(12014h in order to be consistent with the comparisons. 

The relation between both profiles is giyen by 
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where zo is the primordial abundance of heayy elements in the 
Sun and [Fe/H] the metallicity. 

The central star is assumed to haye the same mass as our 
Sun and is a solar metallicity s tar, so that [Fe /H1 = 0. Then, 
2 ” = Zo^g where zo - 0.0149 (lLoddersll2003h . We consider a 
low-mass disk with = O.O 3 M 0 in order to g uarantee that 
the fo rmation of giant planets is not performed (iMiguel et al.l 
I2OI Ih . For the characteristic radius we adopt a yalue of 50 au 
which is in agreement with th e different disk’s obseryations stud¬ 
ied by [Andrews et al.l (1201 Ot) . and for the exponent y we adopt 
a yalue of 1.5 according to other plane tary formation works 
(iRaymond et alJlT005l:lMiguel et al.ll201 Ih . The gas surface den¬ 
sity profile of this disk is a factor of 1 - 1.7 higher than that 
proposed by the Minimun Mass Solar Nebula (MMSN) between 
0 and 15 au. It is worth n oting that the simulations deyeloped 
by IRonco & de Elfal (1201 4ll using density profiles with 7=1.5 
show distinctiye results since true water worlds are formed in the 
habitable zon^ (HZ) of the system . 

Finally, we assumed that the protoplanetary disk presents a 
radial compositional gradient. Thus, the water content by mass 
W(r) assumed as a function of the radial distance r is described 

' The definition of the habitable zone is postponed to Sec. 5 but we 
anticipate that we consider a conse rvative and an optimistic habitable 
zone as in lKoDParapu et af](l201.3al f3i. 
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as: embryos and planetesimals contain 0.001% water inside 2 an, 
contain 0.1% water between 2 an and 2.5 an, contain 5% wa¬ 
ter between 2.5 an and 2.7 an, and 50 % water past 2.7 an. This 
water distribution is the same u sed by iRonco & de Elfal (1201 4ll 
and is based on that proposed bv iRavmond et alJ (l2004l) . Such a 
distribution is assigned to each body in the N-body simulations, 
based on its starting location. Our model does not consider water 
loss during impacts, thus, the final water contents of the resulting 
planets will represent upper limits. 


3. Initial conditions 

As mentioned in Sec. 1, many works that investigate the process 
of planetary formation consider ad hoc initial conditions, partic¬ 
ularly dividing the solid mass of the disk equally between the 
planetesimal and embryo component s. Moreover, some works 
(IChamberil2Q0lt lO’Brien et al.ll2006t) assume that the individ¬ 
ual mass of the planetary embryos has a constant value through¬ 
out the portion of the disk under consideration. However, it is 
worth emphasizing that these conditions represent a particular 
state of the disk. To consider that the solid mass in the inner 
region of the disk is equally distributed between embryos and 
planetesimals could be a reasonable assumption because in fact, 
this depends on the evolution of the system. However, consider¬ 
ing that the population of planetesimals is homogeneously dis¬ 
tributed throughout the inner region of the disk could not be. 
Planetesimal accretion rates are proportional to the surface den¬ 
sity of solids and to the keplerian frequency. These two quantities 
increase closer to the star. Thus, it is expected that inner embryos 
grow quickly accreting all the available planetesimals. Due to the 
fact that embryos can only accrete planetesimals from their feed¬ 
ing zones, which are very narrow as we get closer to the central 
star, we expected a considerable number of small embryos, and 
practically no planetesimals in the inner region at the moment of 
the dissipation of the disk. 

Another consideration that could be not realistic is that all the 
embryos have the same mass in the inner region of the disk. Ba¬ 
sically, the final mass of an embryo (during the gaseous phase) is 
a function of the distance to the star and of the time-scale dissi¬ 
pation of the disk. Brunini & Benvenuto (2008) showed that this 
functional dependence is maximized near the snow line. Thus, it 
is expected that the most massive embryos of the distribution are 
located near the snow line. 

In order to test the sensitivity of the results to the initial dis¬ 
tributions of planetary embryos and planetesimals, we construct 
two different sets of simulations from the same protoplanetary 
disk. The first one uses ad hoc initial conditions for planetesi¬ 
mals and embryos like in iRonco & de Elf^ (l2014t) while the sec¬ 
ond set assumes initial conditions derived from a semi-analytical 
model, which analyzes in detail the evolution of embryos and 
planetesimals during the gaseous phase. 

Since our main goal is to analyze the formation process of 
terrestrial planets and to focus in those of the HZ, our region of 
study extends between 0.5 au and 5 au. Our model assumes an 
inner zone in the disk between 0.5 au and the snow line, and an 
outer zone, which extends between the snow line and 5 au. 

In the following, we describe the generation of initial condi¬ 
tions for both sets of numerical simulations. 

3.1. First set: ad hoc initial conditions 

In this first set of simulations, we assu med the ad hoc init i al con - 
ditions used in the work developed bv iRonco & de Elf^ (l20I4ll . 
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Fig. 1. Initial and final distributions of embryos and planetesimals gen¬ 
erated by the semi-analytical model. Figures a) and b) represent the 
initial and final mass distribution of embryos. At the beginning (a) there 
were 223 embryos between 0.5 au and 5 au, and after 3 Myr of evolution 
(b) there only remain 57 embryos. Figures c) and d) represent the initial 
and final density profile of planetesimals. Figures b) and d) then rep¬ 
resent the final results obtained with the semi-analytical model after 3 
Myr of evolution and the initial conditions for the N-body simulations. 
The vertical lines in grey represent the position of the snow line at 2.7 
au. Color figure is only available in the electronic version. 


Considering a protoplanetary disk of 0.03 Mq and an exponent y 
= 1.5 associated with the surface density profile, the total mass of 
solids in the region of study is 13.66 M^. In agreement with var - 
ious planetary accretion studies such as iKokubo & Idal (Il998l) . 
this mass of solids is divided equally between the planetesimal 
and embryo components. On the one hand, 45 planetary embryos 
are distributed in the study region, 35 of which are in the inner 
zone with masses of 0.06 Mq, and 10 in the outer zone with 
masses of 0.47 Mq. On the other hand, 1000 planetesimals are 
distributed in the disk with masses of 2.68 x 10“^ M® and 0.021 
Mq in the inner and outer zone, respectively. 


3.2. Second set: initial conditions from a semi-analytical 
model 

In this second set of simulations, we generate initial condi¬ 
tions for embryos and plan etesimals using a sern i -analy tical 
model based on the works of iBrunini & Benvenutol (l2008ll and 
iGuilera et alJ (1201 Ol) with the inclusion of some minor improve¬ 
ments. This model analyze the evolution of a protoplanetary disk 
during the gaseous phase and thus, it is possible to obtain more 
realistic initial conditions to be used in a N-body code. A de¬ 
tailed description of this model can be found in Ide Elfa et al.l 
(1201 3h Sec.3. 

In this semi-analytical model, the embryos are separated by 
10 mutual Hill radii and the mass of each embryo corresponds 
to the transition mas s between runaway and oligarchic growth 
(llda & Makinol[T99^ . At the beginning there are 223 embryos 
in the study region, between 0.5 au and 5 au, which present a 
density of 3 gr cm^^ and the planetesimals present a density of 
1.5 gr cm“^ and have radii of 10 km. The population of plan¬ 
etesimals is represented by the density profile of planetesimals. 
In this model, the initial planetesimal size distribution evolves in 
time through planetesimal migration and accretion by embryos. 
Embryos then grow by accretion of planetesimals and also by 
collisions between them. We do not include type I migration 
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in our model and the main reason for neglecting this effect is 
that many quantitative aspects of the type I migration are still 
uncertain as we have already mentioned in Sec. 1. Besides, al¬ 
though planetesimals are affected by the gas drag, as they are 
so big (10 km), this migration is not important to produce sig¬ 
nificant changes in the final amounts of water on embryos and 
in the distribution of water along the planetesimal population. 
We consider that when the distance between two embryos be¬ 
comes smaller than 3.5 mutual Hill radii, they merge into one 
object. We assume perfect inelastic collisions. For mergers be¬ 
tween embryos, we neglected the presence of embryo gaseous 
envelopes because the masses of the embryos are small enough 
not to produce major differences. 

The final configuration obtained after the gas is completely 
dissipated in 3 Myr is also in agreement with that used ad hoc 
in several works as mentioned in Sec. 3.1: we found the same 
proportion for both populations, half the mass in embryos and 
half the mass in planetesimals. Figure [T] shows the results of our 
semi-analytical model for the distribution of embryos and plan¬ 
etesimals for a disk of O.O3M0 at the end of the gas phase. In 
particular, the top left panel represent the mass distribution of 
planetary embryos as a function of the distance from the central 
star used to initiate the evolution of the gas phase with the semi- 
analytical model, while the left bottom panel shows the surface 
density profile of planetesimals used to do the same. The top 
right panel represents the final distribution of embryos obtained 
with the semi-analytical model after 3 Myr of evolution until the 
gas is already gone and the bottom right panel represents the final 
surface density profile of planetesimals. These (right panels) fi¬ 
nal distributions of embryos and planetesimals represent the ini¬ 
tial conditions to be used in our N-body simulations. In this case, 
57 planetary embryos are distributed between 0.5 au and 5 au and 
1000 planetesimals with individual masses of 7.24 x lO^^Me are 
generated using the distribution shown in Fig. Id). 

4. N-body simulations: characteristics 

Our N-body simulations begin when the gas of the disk has al¬ 
ready dissipated. The numerical code used in our N-body simu - 
lations is the MERCURY code developed bv IChambersI (Il999ll . 
We particularly adopted the hybrid integrator, which uses a 
second-order mixed variable simplectic algorithm to treat the in¬ 
teraction between objects with separations greater than 3 Hill 
radii, and a Burlisch-Stoer method for resolving closer encoun¬ 
ters. 

In both sets of simulations, all collisions were treated as 
inelastic mergers that conserve mass and water content. Given 
that the N-body simulations are very costly numerically, we de¬ 
cided to reduce the CPU time considering only gravitational in¬ 
teract ions between embryos and planetesimals jRavmond et alJ 
|2006|) . Thus, planetesimals do not gravitationally interact and 
neither collide with each other. These N-body high-resolution 
simulations allowed us to describe in detail the dynamical pro¬ 
cesses involved during the formation and post evolution stages. 

Since terrestrial planets in our S olar System might 
have formed in 100 Myr - 200 Myr dToubou^l et al.l l20()7t 
[Dauphas & Pourmandll201 lb . we integrated each simulation for 
200 Myr. To compute the inner orbit with enough precision we 
used a 6 day time step. Moreover, the Sun’s radius was artifi¬ 
cially increased to 0.1 au in order to avoid numerical error for 
small-perihelion orbits. 

The orbital parameters, such as initial eccentricities and in¬ 
clinations, were taken randomly considering values lower than 
0.02 and 0.5°, respectively, both for embryos and planetesimals. 



Fig. 2. a) Distributions of embryos used to start the N-body simula¬ 
tions. The squares represent the distribution of embryos obtained with 
the first set of simulations and the circles represent the final results ob¬ 
tained with the semi-analytical model, which represent the second set 
of simulations, b) Surface density profiles used to distribute 1000 plan¬ 
etesimals to start the N-body simulations. The dashed line represents 
the surface density for the first set and the solid line represents the sur¬ 
face density for the second set obtained with the semi-analytical model. 
Color figure is only available in the electronic version. 

In the same way, we adopted random values for the argument 
of pericenter o), longitude of ascending node Q. and the mean 
anomaly M between 0° and 360°. 

Due to the stochastic nature of the accretion process, we per¬ 
formed three simulations SI, S2 and S3 with different random 
number seeds for every set of initial conditions. It is worth men¬ 
tioning that each simulation conserved energy to at least one part 
in lOl 

5. Results 

5.1. General comparative analysis 

The initial distributions associated to planetary embryos and 
planetesimals for both sets of simulations used as initial con¬ 
ditions for the N-body simulations are shown in Figure |3 The 
blue distribution s corres pond to the simulations performed by 
iRonco & de Ella! d2014l) and the red ones correspond to the 
new set of simulations performed from a semi-analytical model. 
Panel a) presents the initial distributions of embryos correspond¬ 
ing to the end of the gas phase and panel b) presents the surface 
density profiles used to distribute 1000 planetesimals in the study 
region, also when the gas is already dissipated. 

As we can see, the embryos inside the snow line of the sec¬ 
ond set of simulations present masses that range from 1 to 2 
times the masses of the embryos inside the snow line of the 
first set. Something similar happens with the embryos beyond 
the snow line. In fact, the masses of the outermost embryos of 
the first set of simulations are at most 5 times greater than those 
of the second set. Another difference evidenced in panel b) is 
that the surface density of planetesimals of both sets are quite 
different in the region inside the snow line since the surface den- 
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Fig. 5. Final configuration of the first and second set of simulations. 
The color scale represent the water content and the shaded regions, the 
optimistic and the conservative HZ. The eccentricity of each planet is 
shown over it, by its radial movement over an orbit. Color figure only 
available in the electronic version. 
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sity profile of planetesimals for the second set of simulations 
is almost zero between 0.5 au and 2.7 au. This is because the 
planetesimals in this region were accreted by the embryos dur¬ 
ing the gaseous phase. Could these differences lead to significant 
changes in the final results of our simulations or not? A discus¬ 
sion on this issue will be addressed in Section 6 . 

Figures[3]and|4]show six snapshots in time on the semimajor 
axis eccentricity plane of the evolution of the S2 and S3 sim¬ 
ulations of the first and the second set, respectively. In general 
terms, the overall progression of both sets of simulations can 
be described as follows. From the beginning, the planetary em¬ 
bryos are quickly excited by their own mutual gravitational per¬ 
turbations. At the same time, the planetesimals, which are not 
self-interacting bodies, significantly increase their eccentricities 
due to perturbations from embryos. In time, the eccentricities of 
embryos and planetesimals increase until their orbits cross and 
accretion collisions occur. Then, planetary embryos grow by ac¬ 
cretion of other embryos and planetesimals and the total number 
of bodies decreases. Finally, at the end of the simulations, the 
entire study region from 0.5 au and 5 au contains between four 
and seven planets with separations ranging from 18 to 53 mutual 
hill radii and a total mass that represents between 55% and 56% 
of the initial mass of solids in the disk which is I 3 . 66 M 0 . 

Between the resulting planets of each simulation, it is possi¬ 
ble to distinguish three classes of special interest; 1 ) the inner¬ 
most planet of the system; 2) the planet (or planets) in the HZ; 
and 3) the most massive planet of the system. Tables [1] and 
show the general characteristics of these three distinctive plan¬ 
ets resulting from SI, S2, and S3 simulations for the first and the 
second set of initial conditions, respectively, and Fig.|5]shows the 
final configuration of all the simulations performed. The main 
similarities and differences obtained in these simulations for the 
resulting planets as well as the remaining planetesimal popula¬ 
tion are discussed throughout the following subsections. 


5.2. Removal of embryos and planetesimals 

It is worth mentioning that while the total mass of the planetes¬ 
imal population starting the N-body simulations is almost the 
same for both sets of initial conditions (~ 6.9M@ and 7.2M® for 
the first and the second set, respectively), the initial surface den¬ 
sity profile of planetesimals is very different in both cases (see 
Fig.|2]r). In fact, for the first set of simulations, 31% of the mass 
of the planetesimal population is located inside the snow line, 
while the remaining 69% is distributed beyond 2.7 au. On the 
other hand, for the second set of initial conditions, only 10% of 
the mass of the planetesimal population is distributed inside the 
snow line, while the remaining 90% is located beyond 2.7 au. 
Thus, some differences should be observed in the evolution of 
the planetesimal population between both sets of initial condi¬ 
tions. 

After 200 Myr of evolution, many embryos and planetesi¬ 
mals were removed from the disks of both sets of simulations. 
In the S2 simulation of the first set of initial conditions, the mass 
of planetary embryos and planetesimals that still remain in the 
study region, between 0.5 au and 5 au, is ~ 7M^ and 0.43M®, 
respectively. We consider that this remaining planetesimal mass, 
which represents 6.3% of the initial planetesimal mass, is not 
significant to produce important changes in the orbital and phys¬ 
ical characteristics of the resulting planets. Similar results were 
found in the other simulations of this set. 

In the S3 simulation of the second set of initial conditions, 
the mass that still survive after 200 Myr is 7.49M^ in embryos 
and 1.08M® in planetesimals. Thus, the remaining mass of plan¬ 
etary embryos in the study region is comparable in both sets 
of initial conditions. However, the mass of remaining planetes¬ 
imals associated to the second set of simulations, which repre¬ 
sents 15% of the initial planetesimal mass, is a factor of ~ 2.4 
greater than that obtained from the first set of initial conditions. 
This situation is repeated in the rest of the simulations of this set. 

From a quantitative point of view, the populations of em¬ 
bryos at the end of a simulation obtained from both sets of initial 
conditions result to be comparable, while the remaining popu¬ 
lations of planetesimals show the main discrepancies. In order 
to understand such differences, it is very important to analyze 
the main mechanisms of removal of planetesimals in our simu¬ 
lations. On the one hand, in the S2 simulation of the first set, 
3 . 78 M 0 in planetesimals, which represents 55% of the initial 
mass in planetesimals, collide with the central star, remains out¬ 
side the study region (beyond 5 au) or is ejected from the disk. 
On the other hand, in the S3 simulation of the second set 3.83M0, 
which represents 53% of the initial mass in planetesimals, col¬ 
lide with the central star, remains outside the study region or is 
ejected from the disk. Then, the removed mass in planetesimals 
is also comparable in both sets of simulations. However, the ac¬ 
creted mass in planetesimals (by embryos) represents 37.9% of 
the initial mass in planetesimals for the S2 simulation of the first 
set and represents 32.3% of the initial mass in planetesimals for 
the S3 simulation of the second set. Similar results were found 
for the rest of the simulations. While the removed mass in plan¬ 
etesimals is comparable in both sets and range between ~ 53% 
and 55% of the initial planetesimal mass for the second and 
the first set, the percentage of the accreted mass in planetesi¬ 
mals range between 37.9% and 41% for the first set and between 
31.7% and 32.3% for the second set. This difference in the ac¬ 
creted mass in planetesimals between both sets of simulations 
reveals distinct accretion histories of the final planets. The ac¬ 
cretion history of the planets is of great interest because the final 
characteristics and properties of the planets depend on it. Thus, 
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Fig. 3. Evolution in time of the S2 simulation of the first set. The blue and light-blue shaded areas represent the conservative and the optimistic HZ 
respectively and the blue and light-blue curves represent curves of constant perihelion and aphelion, both for the conservative and the optimistic 
HZ. Planetary embryos are plotted as colored circles and planetesimals are plotted with black dots. The color scale represents the fraction of water 
of the embryos respect to their total masses. In this case there are two planets within the optimistic HZ with masses of I.I 9 M 0 and I. 65 M 0 , 
respectively. They present 4.51% and 39.48% of water content by mass, which represents 192 and 2326 Earth’s oceans, respectively. Color figure 
is only available in the electronic version. 


Table 1. General characteristics of the most distinctive planets formed in the SI, S2, and S3 simulations of the first set, using ad hoc initial 
conditions, is the initial and af is the final semimajor axis. The initial semimajor axis of a planet represents the initial location of the most 
massive embryo from which the planet formed. < e >iooMyr and < ' >iooMyr represent the mean eccentricity and inclination, respectively, in the 
last 100 Myr of evolution. M{M^) represents the final mass, VP(%) the final percentage of water by mass and Tlgi the timescale in Myr associated 
with the last giant impact produced by an embryo. 


Simulation 

Planet 

ai(au) 

flf(au) 

< e >100Myr 

< i >100Myr 

M(Me) 

W(%) 

Tlgi (Myr) 

SI 

Innermost 

0.68 

0.62 

0.07 

4° .25 

1.57 

5.47 

56.1 


Planet a in HZ 

4.52 

1.41 

0.06 

2°.39 

2.21 

32.55 

35 


Most massive 

4.52 

1.41 

0.06 

2°.39 

2.21 

32.55 

35 

S2 

Innermost 

0.78 

0.47 

0.14 

5°.86 

0.6 

7.12 

5.55 


Planet b in HZ 

1.30 

0.90 

0.08 

3°.56 

1.19 

4.51 

35 


Planet c in HZ 

3.13 

1.35 

0.06 

2°.86 

1.65 

39.48 

22 


Most massive 

3.25 

3.35 

0.05 

2°.20 

2.15 

43.66 

96.8 

S3 

Innermost 

0.50 

0.63 

0.08 

4° .04 

1.13 

6.66 

11.45 


Planet d in HZ 

0.64 

0.98 

0.07 

3°.64 

0.66 

8.10 

6.00 


Planet e in HZ 

3.43 

1.53 

0.05 

2°.23 

2.27 

32.46 

71.9 


Most massive 

4.69 

2.77 

0.04 

1°.47 

3.08 

47.86 

44.3 


initial conditions obtained from more realistic models allow us 5.3. Dynamical friction 
to obtain more reliable physical properties of the resulting plan¬ 
ets like, for example, a more reliable final water distribution. From the beginning, the dynamical friction seems to be a signif¬ 
icant phenomenon in both sets of simulations. This dissipative 
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Fig. 4. Evolution in time of the S3 simulation of the second set. The blue and light-blue shaded areas represent the conservative and the optimistic 
HZ respectively and the blue and light-blue curves represent curves of constant perihelion and aphelion, both for the conservative and the optimistic 
HZ. Planetary embryos are plotted as colored circles and planetesimals are plotted with black dots. The color scale represents the fraction of water 
of the embryos respect to their total masses. Two planets remain within the optimistic HZ with masses of 1.31 and l.hSM®, respectively. They 
present 7.49% and 24.38% of water content by mass, which represents 366 and 1437 Earth’s oceans, respectively. Color figure is only available in 
the electronic version. 


Table 2. General characteristics of the most distinctive planets formed in the SI, S2, and S3 simulations of the second set, using initial conditions 
from a semi-analytical model, is the initial and Uf is the final semimajor axis. The initial semimajor axis of a planet represents the initial location 
of the most massive embryo from which the planet formed. < e >iooMyr and < i >iooMyr represent the mean eccentricity and inclination, respectively, 
in the last 100 Myr of evolution. M(M^) represents the final mass, 'W{%) the final percentage of water by mass and Flgi the timescale in Myr 
associated with the last giant impact produced by an embryo. 


Simulation 

Planet 

a/au) 

flf(au) 

< e >100Myr 

< i >100Myr 

M(Me) 

W{%) 

Tlgi (Myr) 

SI 

Innermost 

0.54 

0.44 

0.12 

8°.76 

0.45 

13.03 

74.24 


Planet a in HZ 

1.05 

0.78 

0.05 

2°.36 

2.21 

8.75 

51.6 


Planet b in HZ 

1.49 

1.38 

0.06 

3°.02 

1.18 

24.30 

28.2 


Most massive 

1.05 

0.78 

0.05 

2°.36 

2.21 

8.75 

51.6 

S2 

Innermost 

0.87 

0.56 

0.07 

3°.93 

1.33 

10.93 

198.8 


Planet c in HZ 

0.73 

1.04 

0.04 

2 ° .01 

2.00 

10.76 

20.5 


Planet d in HZ 

2.73 

1.70 

0.05 

2° .09 

1.67 

33.63 

90.5 


Most massive 

0.73 

1.04 

0.04 

2°.01 

2.00 

10.76 

20.5 

S3 

Innermost 

0.87 

0.58 

0.06 

2°.47 

1.18 

9.90 

34.8 


Planet e in HZ 

1.33 

0.88 

0.06 

2° .09 

1.37 

7.49 

18.3 


Planet/ in HZ 

1.14 

1.37 

0.05 

1°.67 

1.65 

24.38 

59.7 


Most massive 

2.92 

2.55 

0.03 

1°.55 

2.36 

39.74 

62.7 


force dampens the eccentricities and inclinations of the large 
bodies embedded in a swarm of smaller bodies represented by 
planetesimals. Figure |6] shows the evolution of the eccentricities 


and inclinations of the innermost planet and the most massive 
planet in S3 of the first set and in S3 of the second set. The 
masses of the innermost planets are comparable and the same 
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happens with the most massive planets (see Tables [T] and |2]i. 
Thus, we find similar results for both of them. 

On the one hand, the innermost planets of 1.1 3 M 0 and 
1 . 18 M 0 show mean values of the eccentricity of 0.08 and 0.06 
for the first and the second set, respectively, while the mean val¬ 
ues for the inclinations are 4.04° and 2.47°. On the other hand, 
the most massive planets of 3 .O 8 M 0 and 236M^ present mean 
values of the eccentricity of 0.04 and 0.03 for the first and the 
second set, respectively, while the mean values for the inclina¬ 
tions are 1.47° and 1.55°. This means that, for both kinds of sim¬ 
ulations, the effects of dynamical friction still prevail over larger 
bodies. It is worth noting that both planets are initially located in 
regions embedded in planetesimals, therefore, it is expected to 
see the same effects in both of them. 

Similar results are obtained with the comparison between S2 
of the first set and SI of the second set. In this particular sce¬ 
nario, it is worth noting that the most massive planet in S1 of the 
second set is not located in a region embedded in planetesimals 
at the beginning of the integration. However, this planet has of 
the order of 5 x 10^ close encounters with planetesimals during 
its evolution, which allow to keep the low values associated to its 
eccentricity and inclination. Therefore, the fact that the planet is 
located at the beginning in a region without planetesimals does 
not implies that the effect of dynamic friction is absent through¬ 
out its evolution. 

An interesting scenario to analyze is represented for the S2 
simulation of the second set. In this particular case, the inner¬ 
most and the most massive planets started the simulation in a 
region of the disk which is not embedded in planetesimals at the 
beginning of the integration. It is important to note that the final 
masses of such planets have comparable values. However, the 
mean values of the eccentricity and inclination of the innermost 
planet are a factor of ~ 2 greater than those associated to the 
most massive planet. Such a difference is due to the fact that the 
most massive planet suffered one order of magnitude more close 
encounters with planetesimals that the innermost planet through¬ 
out the evolution. Thus, the results of our simulations suggest 
that the dynamical friction is an inefficient mechanism for the 
innermost planets. 


5.4. Planets in the habitable zone 


The formation of terrestrial planets in the HZ of the system and 
the analysis of their water contents are topics of special inter¬ 
est in our work. The HZ is defined as the range of heliocen¬ 
tric dist ances at which a pl anet can retain liquid water on its 
surface (i Kasting et al. 19931) . For the case of the Solar System, 
iKopparapu et al.l2013allbh proposed that an optimistic HZ is be¬ 
tween 0.75 au and 1.77 au and a conservative HZ is between 
0.99 au and 1.67 au. However, it is worth noting that the fact 
that a planet is in the HZ does not guarantee that there may be 
developed life. Planets with very eccentric orbits may pass most 
of their periods outside the HZ, not allowing long times of per¬ 
manent liquid water on their surfaces. To avoid this problem, we 
considered that a planet is in the optimistic HZ if it has a per¬ 
ihelion q > 0.75 au and a aphelion Q < 1.77 au and is in the 
conservative HZ if it has a perihelion q > 0.99 au and a aphelion 
Q < 1.67 au. In this work, we considered it is sufficient that a 
planet is in the optimistic HZ to be potentially habitable. On the 
other hand, we considered the water contents to be significant 
when they are similar to that of the Earth. The mass of water on 
the Earth surface is 2.8 x 10 M^, which is defined as 1 Earth 

ocean. Taking i nto account the water content in the Ea rth mantle, 
several studies dLecuver & Gilletlll998l : lMartvll2012h suggested 
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Fig. 6. Evolutions of eccentricities, a), b), and inclinations, c), d), as 
a function of time for the innermost planet (black curve) and the most 
massive planet (red curve) of the S3 simulations and S3 simulations 
in the first and in the second set, respectively. The dynamical friction 
effects prevail over the larger bodies which evolve, in both sets of sim¬ 
ulations, in regions of the disk embebbed in a swarm of planetesimals. 
Color figure is only available in the electronic version. 



that the present-day water content on Earth should be ~ 0.1% 
to 0.2% by mass, which represents between ~ 3.6 to 7.1 Earth 
oceans. Moreover, some works suggest that the primitive man¬ 
tle of t he Earth could have contained between 10 and 50 Earth 
oceans dAbe et al.l2000t[Righter & Drake! 19^^ . although this is 
still under debate. 

We developed three different N-body simulations for every 
of the two sets of initial conditions. It is worth noting that all of 
them form planets within the optimistic HZ. Eor the first set of 
initial conditions, our simulations produce a total of five plan¬ 
ets in the HZ with masses from 0.66 M 0 to 2.27 M 0 and final 
water contents ranging from 4.5% to 39.48% of the total mass, 
which represents from 192 to 2326 Earth oceans. Eor the second 
set, our simulations form a total number of six planets in the HZ 
with masses ranging from 1.18 M 0 to 2.21 M 0 and final water 
contents from 7.5% to 33.63% of the total mass, which repre¬ 
sents from 427 to 2006 Earth oceans. 

Eigure|2] shows the evolution of the semimajor axis, the per¬ 
ihelion and the aphelion of the planets that remain in the HZ at 
the end of the simulation for both sets of initial conditions. As 
we can see, most of the planets are kept within the conservative 
HZ almost throughout the whole integration, and their semima¬ 
jor axes do not change significantly in the last 100 Myr of evolu¬ 
tion. Planet b (S2 simulation) in the first set of initial conditions 
as well as planets a (SI simulation) and e (S3 simulation) in the 
second set remain in the limits of the optimistic HZ, with perihe- 
lions entering and leaving it due to oscillations in the eccentricity 
of their orbits. Since the definition of HZ is not accurate, we will 
consider that these planets are targets of potential interest. Planet 
d (S2 simulation) in the second set of initial conditions reaches 
maximum values of the eccentricity of ~ 0.15. This planet shows 
changes in its aphelion, which carry it in and out of the HZ at the 
end of the simulation. The incident stellar flux on such a planet 
may have considerable variat ions between perihel i on an d aphe¬ 
lion. The study developed by IWilliams & PollardI (l2002h shows 
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Fig. 7. Evolution in time of the semimajor axis (black), the perihelion (red) and the aphelion (orange) of the planets that remain within the HZ 
of each simulation in both, the first (left panel) and the second set (right panel) of simulations. The blue and light-blue shaded areas represent the 
conservative and the optimistic HZ respectively and the dashed grey line represents the position of the snow line. Color figure is only available in 
the electronic version. 


that, provided that an ocean is present to act as a heat capacitor, 
it is primarily the time-averaged flux that affects the habitabil¬ 
ity over an eccentric orbit. Given that planets on eccentric orbits 
have higher average orbital flux, planet d (S2 simulation) in the 
second set may maintain habitable conditions. 

5.4.1. Highly water-rich planets 

As the reader can see in Tables 1 and 2, simulations developed 
from both sets of initial conditions form planets in the HZ with 
a wide range of final water contents. In particular, some of them 
show very high water contents. In fact, planets a, c and e in the 
first set of simulations show water contents of 32.6%, 39.5% 
and, 32.5% by mass, respectively, while planets b, d and / in 
the second set present 24.3%, 33.6%, and 24.4% water by mass, 
respectively. It is worth noting that the evolution of such water- 
rich planets has different characteristics for both sets of simula¬ 
tions. Tables 1 and 2, and Figure|7]give us an interesting vision 
concerning this point. In fact. Table 1 and Figure |7] (left panel) 
show us that potentially habitable planets a, c and e of the first 
set of simulations come from beyond the snow line. The migra¬ 
tion suffered by these planets is due to the strong gravitational 
interaction between planetary embryos and planetesimals in the 
outer disk. Given their initial locations, it is worth noting that 
an important percentage (31%-36%) of the very high final water 
contents of such planets results to be primordial. Moreover, our 
simulations indicate that the remaining water content of those 
planets is acquired by impacts of planetesimals and planetary 


embryos in equal parts. Thus, taking into account the primordial 
and acquired water contents, we infer that the embryos play the 
most important role to provide the water to the potentially habit¬ 
able planets a, c and e of the first set of simulations. 

Table 2 and Figure |7] (right panel) show us that potentially 
habitable planets b, d and / of the second set of simulations 
present different characteristics. On the one hand, planet d starts 
the simulation beyond the snow line, like the potentially habit¬ 
able planets a, c and e of the first set. In the same way, 41.4% of 
the high final water content of this planet is primordial, being a 
consequence of its initial location. Moreover, the remaining wa¬ 
ter content is provided equally by impacts of embryos and plan¬ 
etesimals associated to the outer disk. Thus, taking into account 
the primordial water content, we find that the embryos are the 
main responsible to provide the water to planet d of the second 
set of simulations. On the other hand, planets b and / present a 
semimajor axis which do not change significantly from the be¬ 
ginning of the integration. This means that such planets grow 
near their initial positions and do not suffer significant migra¬ 
tion. Given their initial locations, the primordial water content 
of the planets b and / is negligible. In fact, for both planets, the 
final water content is fully acquired throughout evolution by im¬ 
pacts of planetesimals and planetary embryos associated to the 
outer disk. This can be seen in the bottom panel of Fig. [8]where 
the feeding zones of the planets are represented. In particular, 
planets b and / accreted two and one embryos from beyond the 
snow line respectively, which significantly increased their final 
water contents. However, planetesimals provide 53% and 44% of 
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the final water contents of planet b and/. Thus, for these cases, 
both embryos and planetesimals are equally responsible to pro¬ 
vide the final water contents of such planets. 

Although these highly water-rich planets remain within the 
optimistic h abitable zone, their potential habitability is still un- 
der debate (lAbbot et alJ 120121: lAliberll 1201 4t iKitzmann et alJ 
I2OI5I) and that analysis is beyond the scope of this work. How¬ 
ever we consider water worlds as a particular and very interesting 
kind of exoplanets. 

5.4.2. Earth-like planets 

Our simulations also form Earth-like or terrestrial-type planets 
within the optimistic HZ. We consider the Earth is comparable 
to these planets not only because they present masses similar to 
that of the Earth, but rather because of the features that have to do 
with the dynamics of their formation. These Earth-like planets in 
our simulations form in situ. There is no significant migration of 
the accretion seeds that end up forming such planets. Moreover, 
their feeding zones, unlike what happens with highly water-rich 
planets, are restricted to the inner region of the disk (see Eig. 
EJ. This kind of exoplanets results to be the most interesting one 
from an astrobiological point of view since, as they present sim¬ 
ilarities in terms of the dynamics of the formation of the Earth, 
they are more likely to develop life. It is worth mentioning that 
these Earth-like planets formed in our simulations present sev¬ 
eral times the Earth’s water content. This is directly related to 
two important facts: on the one hand, as we mentioned before, 
our model does not consider water loss during impacts, thus, the 
final water contents represent upper limits. On the other hand, 
our formation scenario without gas giants is particularly con¬ 
ducive to a more efficient water delivery from external planetesi¬ 
mals. Giant planets could act as a barrier preventing that many of 
the external planetesimals, which are now accre ted by the plan¬ 
ets in the HZ, reach the inner zone of the disk (iRavmond et alJ 
120091) . 

Two of the five planets formed within the optimistic HZ in 
the first set and three of the six planets that remain in the HZ 
of the second set are Earth-like planets. The Earth-like planets in 
the first set (see Table 1: planet b in S2 and planet d in S3) present 
masses of 1.1 9 M 0 and O. 66 M 0 and water contents of 4.51% and 
8.10% by mass, respectively, which represent 192 and 190 Earth 
oceans. The Earth-like planets in the second set (see Table 2: 
planet a in SI, planet c in S2 andplanet e in S3) present masses of 
2.2 IM 0 , 2 M 0 and 1 . 37 M 0 , and water contents of 8.75%, 10.76% 
and 7.49% by mass, respectively, which represent 692, 768 and 
366 Earth oceans. 

The evolution of these Earth-like planets presents similar 
characteristics in both sets of simulations unlike what happens 
with the water worlds already described. Eigure [T] shows that 
planets b and d in the first set, and planets a, c and e in the second 
set, start their formation in the inner zone of the disk, particularly 
inside 1.5 au, and they evolve without exceeding this limit. This 
means that the percentage of water that they have achieved dur¬ 
ing their evolution is not primordial. In fact, a detailed analysis of 
their collisional histories reveals that planetesimals are the only 
population responsible for their final amounts of water. More¬ 
over, none of them accreted embryos from beyond the snowline. 
However there is a difference regarding what class of planetes¬ 
imals these planets have accreted. Erom now on we will call 
planetesimals class-D to those planetesimals which were orig¬ 
inally located inside the snow line and are dry, and we will call 
planetesimals class-W to those planetesimals which were origi¬ 
nally located beyond the snow line and present half of their mass 


in water. Planets b and d of the first set accreted 185 and 100 
planetesimals of the total of planetesimals, of which only 5 were 
class-W. This represents only 2.7% and 5% of the total amount 
of accreted planetesimals. In contrast, planets a, c and e of the 
second set accreted 65, 69 and 41 planetesimals of the total ac¬ 
creted planetesimals, of which 53, 59 and 28 were class-W plan¬ 
etesimals. This represents 81.5%, 85.5% and 68 % of the total 
amount of accreted planetesimals for planets a, c and e, respec¬ 
tively. Thus, the class-W planetesimals accreted for both sets of 
simulations is quite difference. Nevertheless it is necessary to 
remember that the masses of the class-W planetesimals are dif¬ 
ferent in both sets of simulations. We can note that the mass of 
water accreted by the Earth-like planets considering only class- 
W planetesimals in the second set is between 1.9 and 4 times the 
mass of water accreted by the planets of the first set. 

An important topic to consider is the comparison between the 
timescale for an Earth-like planet to accrete half of the class-W 
population of planetesimals (fso.dass-w), the timescale to accrete 
half of the class-D planetesimals (fso.ciass-o) and the timescale for 
the same planet to achieve at least 50% of its final mass 
Table [3] shows the main characteristics of the Earth-like planets 
formed within the HZ in both sets of simulations and shows that, 
for all of them, fso.ciass-w is greater than fso.dass-D, greater than 
7 ’ 5 o% and even more, is greater than T75%- This means that class- 
W population of planetesimals was accreted significantly later 
than class-D planetesimals during the formation process and af¬ 
ter the planet achieved half of its final mass. Moreover, class-W 
planetesimals were accreted after the planet achieved 75% of its 
final mass. This late arrival of the class-W planetesimals dur¬ 
ing the formation of terrestrial planets may favor water reten- 
tion on the surface of the planet brought by these planetesimals. 
lO’Brien et al.l (1201 4l) studied the water delivery and giant im¬ 
pacts on_ten;estrial_planets in the “Grand Tack” scenario. Particu¬ 
larly, [oiBrien|et^ d2014l) analyzed the wate r delivery of primi- 
tive planetesimals (“C-Type” planetesimals in i Walsh et al.lbolU) 
which were initially located between the orbits of Jupiter and 
Saturn, or beyond. Although the Grand Tack scenario is quite 
different from the one we proposed since our planetary systems 
do not include giant planets, we find similar results concerning 
timescales of primitive planetesimals accretion. 

Moreover, we find that the Earth-like planets of the first and 
the second set are still accreting water due to class-W planetesi¬ 
mals after ~ 90 Myr of evolution. Although planets b and d of the 
first set only accreted 5 class-W planetesimals, 3 of them were 
accreted after 90 Myr, particularly between 110 and 170 Myr. 
These planetesimals provided 58.7% and 39.3% of the total ac¬ 
creted water by planets b and d, respectively, which represent 
2.64% and 3.2% of their final masses. Eor the second set, planet 
a accreted 10 class-W planetesimals between 93 and 191 Myr 
which provided 18.7% of the total accreted water, and represent 
1.64% of its final mass. Planet c accreted 20 class-W planetes¬ 
imals between 91 and 199 Myr which provided 33.6% of the 
total accreted water and represent 3.62% of its final mass and 
finally planet e accreted 11 class-W planetesimals between 92 
and 192 Myr which provided 38.8% of the total accreted water 
and represent 2.91% of its final mass. This analysis may be of 
particular interest for the habitabl e zone plan e ts form ed around 
M dwarfs regarding the work of iTian & Idal (l2015l) . Although 
the Earth-mass planets found bv lTian & Idal (120151) in M dwarf 
habitable zones may have lost their water content through dis¬ 
sociation and hydrodynamic scape during the first ~ 90 Myr of 
evolution, they could still accrete ice planetesimals and retain 
significant amounts of water on their surfaces. 
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Table 3. General characteristics of the Earth-like planets formed within the HZ in both sets of simulations. MfM®) represents the final mass, 
tso.ciass-D represents the timescale in Myr for a planet to accrete half of the class-D planetesimals that will accrete, tso.ciass-w represents the timescale 
in Myr for a planet to accrete half of the class-W planetesimals that will accrete, represents the timescale in Myr for a planet to achieve 50% 
of its final mass and Tji% represents the timescale in Myr for a planet to achieve 75% of its final. Planet b of the first set shows the same Tso* and 
Tis%- This means that the planet suddenly accreted enough mass to go from having less than 50% of its final mass to have more than 75% of its 
final mass. The same happens with planet a of the second set. 


Simulations 

Planet 

M(M 0 ) 

W(%) 

Go.ciass-D (Myr) 

Go.ciass-w (Myr) 

Tsofo (Myr) 

Ti 5 % (Myr) 

Eirst Set 

b 

1.19 

4.51 

21.9 

81.3 

35 

35 


d 

0.66 

8.10 

21.8 

613 

2.9 

13.4 

Second Set 

a 

2.21 

8.75 

21.8 

56.6 

51.7 

51.7 


c 

2.00 

10.76 

39.8 

66.9 

6 

20.5 


e 

1.37 

7.49 

27.1 

67.8 

4.4 

18.4 



0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

Semimajor axis (au) 


Fig. 8 . Feeding zones of the planets that remain in the HZ of S1, S2 and 
S3 simulations of the first (top panel) and the second set of simulations 
(bottom panel). The y axis represents the fraction of each planet’s final 
mass after 200 Myr of evolution. Color figure is only available in the 
electronic version. 

6. Discussion 

We found that at the end of the evolution of our simulations, 
five and six planets are produced in the HZ of the system from 
the first and second set of initial conditions, respectively. For the 
first set, planets of the HZ present masses ranging from 0.66M@ 
to 2.27M0, with a median value of 1.65M0. For the second set, 
the mass of the planets formed in the HZ ranges from 1.1 8 M 0 to 
2 . 2 IM 0 , with a median value of I. 66 M 0 . The timescales associ¬ 
ated to the last giant impact received for the planets of the HZ 
are quite similar for both sets of simulations, with mean values 
of 35 Myr and 40 Myr for the first and second set, respectively. 
However, not always the last giant impact means that the planet 
presents the majority of it final mass. Figure |3 shows the mean 
time values that the planets in the HZ need to reach 50%, 75%, 
and 90% of their final mass for both sets. The values are aver¬ 
ages for all planets in the HZ and, as we can see, the results are 
also quite similar for both sets. 

A simple analysis of the feeding zones of the HZ planets 
(Figure [3) reveals that the planets in the first set of simulations 
accrete a higher percentage of material from the outer disk that 
the potentially habitable planets produced in the second set. In 
fact, for the first set, the final water contents range from 4.5% to 
39.5% by mass, with a median value of 32.5%. For the second 



Fig. 9. Timescales for planets in the HZ to reach a given fraction (50%, 
75%, or 90%) of their final masses for the set of simulations with ad hoc 
initial conditions and for the set of simulations with initial conditions 
from a semi-analytical model. These values are averaged for all the HZ 
planets in each simulation. 


set, planets in the HZ have final water contents ranging from 
7.5% to 33.6% by mass, with a median value of 17.5%. 

Our simulations form “Earth-like planets” and “highly 
water-rich planets” or “water worlds” in the HZ. All the Earth¬ 
like planets, in both sets, grow in situ and their feeding zones are 
restricted to the inner disk. Thus, their final water contents are 
not primordial, and a detailed analysis shows that planetesimals 
are the only population responsible for their final amounts of 
water. In addition, an analysis of the accretion timescales shows 
that all these planets accreted half of the class-W planetesimals 
after they have achieved 75% of their final masses. Moreover, 
they continue accreting class-W planetesimals after 90 Myr of 
evolution. On the other side, the Earth-like planets in both sets 
of simulations accreted different ammounts of class-W planetes¬ 
imals. Particularly, the mass of water accreted by the Earth-like 
planets of the second set is greater than the one accreted by the 
first set, by a factor of 1.9 - 4. The highly water-rich planets are 
planets whose accretion seeds come from the outer zone of the 
disk (planets a, c and e of the first set and planet d of the sec¬ 
ond set) or grow near their initial positions but have accreted 
embryos from beyond the snow line. The sourses of water suply 
are also different in both sets. Planetary embryos are the most 
responsible for the final water contents of the water worlds of 
the first set but for two of the three planets of the second set, the 
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final water contents are provided by embryos and planetesimals 
in equal parts. 

However, it is worth mentioning that e mbryos and planetesi- 
mals might not be the only source of water. Ilzidoro et all (l2013h 
suggested that a compound model that considers water accretion 
through collisions of protoplanetary bodies and the contribution 
of the water locally adsorbed from the nebula is more efficient in 
the amount and time of the delivery of water to Earth-like plan¬ 
ets. 

Considering more realistic initial conditions to include in the 
N-body simulations is important to get a more realistic accre¬ 
tion history of the hnal planets. However, it is worth noting that 
the semi-analytical model used to generate these initial condi¬ 
tions contains several simplifications which might affect the re¬ 
sults and could lead to different hnal conhgurations for the dis¬ 
tribution of embryos and planetesimals. Particularly, our model 
neglects the presence of embryo gaseous envelopes because the 
masses of the embryos are small enough not to produce major 
differences. Moreover, the semi-analytical model does not in¬ 
clude the effects of planetesimal fragmentation or pebble accre¬ 
tion. On the other hand, we do not include type 1 migration in our 
model because many aspects of this phenomenon are still under 
debate. Computing the N-body interactions for the embryos dur¬ 
ing the gaseous phase could also change their spacial distribu¬ 
tion. Regarding the planetesimal distribution, different sizes of 
planetesimals could lead to different hnal planetesimal surface 
densities. 

A detailed knowledge about the accretion history of the plan¬ 
ets of a given system is very important in order to analyze sev¬ 
eral physical characteristics. In particular, all our simulations 
are developed assuming perfectly inelastic collisions that con¬ 
serve mass and water content. However, a more realistic treat¬ 
ment of the collisions can lead us to a better understanding 
about the differentiation of a planet, composition, formation of 
structures such as a core and a mantle, as well as to quantify 
its abundance of water (iMarcus et al.l l2010l : IChambera 120131 
iBonsor et al.ll2015ll . At the same way, a better knowledge about 
the accretion history of a planet allows us to und erstand the evo¬ 
lution of its atmosphere dSchlichting et alJ2015h . In fact, impacts 
due to planetary embryos and planetesimals can generate a rel¬ 
evant atmospheric mass loss during the formation and evolution 
process. 

7. Conclusions 

We carried out N-body simulations aimed at studying the plane¬ 
tary formation process and water delivery around Sun-like stars 
in low-mass disks and in absence of gas giants. The main goal 
of our investigation is to study the sensitivity of the results to 
the particular initial distribution addopted for planetesimals and 
planetary embryos aft er the gas dissipation. Th e hrst set of sim¬ 
ulations was based on iRonco & de Elfal (l2014t) and assumes ad 
hoc initial conditions while the second set of simulations derives 
initial conditions from a semi-analytical model, which is capable 
of analyzing the evolution of a planetary system during the gas 
phase. Our main results suggest: 

- The number of planets formed in the HZ of the system is 
not sensitive to the initial conditions proposed in both sets of 
simulations. Moreover, the masses of the planets of the HZ 
do not show signihcant differences in both sets. 

- The main differences observed between both sets of simu¬ 
lations are associated to the accretion history of the planets 
of the HZ. Primarily, such discrepancies are related to the 
accretion of water-rich material. 


- Earth-like planets form in situ and their feeding zones are 
restricted to the inner region of the disk and these Earth-like 
planets present several Earth’s water contents. 

- Highly water-rich planets present different characteristics in 
both sets: these are planets whose accretion seeds come from 
the outer zone of the disk or are planets which grow near their 
initial positions but have accreted embryos from beyond the 
snow line, meaning that their feeding zones extend to the 
outer region of the disk. 

- Our simulations indicate that planetary embryos are the most 
responsible for the final water contents of the water worlds in 
the HZ of the first set, while planetesimals play a secondary 
role. Of the three water worlds formed in the HZ from the 
second set of initial conditions, only one of them shows a 
similar behavior to that previously described. However, the 
primordial water content of two of such planets is negligi¬ 
ble and their hnal water contents are provided by impacts of 
embryos and planetesimals in equal parts. 

- Almost all the simulations formed Earth-like planets within 
the HZ. As they grow in situ, their hnal water contents are 
not primordial, and planetesimals are the only population re¬ 
sponsible for their hnal amounts of water. In addition, all 
these planets accreted half of the class-W planetesimals af¬ 
ter they have achieved 75% of their hnal masses. Moreover, 
they continue accreting class-W planetesimals after 90 Myr 
of evolution. These results show that the water delivery is a 
late process in the evolution of an Earth-like planet and this 
is of great interest because it can facilitate the retention of 
water and other volatiles on the surface. 

Einally we conclude that the inclusion of more realistic ini¬ 
tial conditions in N-body simulations of terrestrial-type planet 
formation is crucial to do a more realistic analysis of the ac¬ 
cretion history of the planets resulting from the formation of a 
planetary system. 
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